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INTRODUCTION 


The  objective  of  this  study  is  to  demonstrate  the  utility  of 
the  phase  retrieval  concept  for  determining  the  wavefront  aberrations 
of  an  optical  system  from  image  intensity  data.  Phase  retrieval 
techniques  use  only  the  focal  plane  detector  array  to  estimate  the 
wave  aberration  function  6  of  an  active  optical  system.  The  esti¬ 
mate  would  then  be  used  to  derive  control  signals  to  align  or 
maintain  alignment  of  the  optical  system.  See  Figure  1. 


object 


adaptive  detector  array  image 

optics 


Figure  1.  Phase  retrieval  in  an  adaptive  imaging  system 

In  a  previous  study  (F30602-77-C-0176)  we  showed  that  phase 
retrieval  is  a  promising  technique  to  maintain  alignment  of  an  active 
optical  system.  However,  most  of  the  study  was  theoretical  in  nature 
and  the  best  results  were  obtained  with  point  sources,  little  or  no 
noise  on  the  detector  output,  and  adequately  sampled  data  (at  or 
higher  than  the  Nyquist  sampling  rate) .  In  this  follow-on  effort 
we  are  charged  to  thoroughly  simulate  the  algorithm's  performance, 
particularly  with  respect  to  noise,  undersampling  (because  of  the 
finite  detector  size),  and  extended  objects.  We  are  also  charged  to 
show  how  phase  retrieval  can  be  developed  into  a  system  level  proof- 
of-concept  demonstration,  where  the  system  is  to  be  specified  by  the 
government. 
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In  this  report  we  review  the  phase  retrieval  concept,  show  the 
effects  of  undersampling  and  noise,  present  algorithms  for  extended 
objects,  show  several  ways  in  which  the  algorithm  can  be  speeded  up, 
introduce  the  concept  of  phase  diversity  imaging  (as  a  research 
by-product  of  our  extended  object  algorithm) ,  give  results  of  an 
experiment  that  attempts  to  verify  the  extended  object  concept, 
and  present  the  results  of  a  series  of  blind  tests  of  the  phase 
retrieval  algorithm  when  the  data  is  noisy  and  undersampled. 
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REVIEW  OF  THE  PHASE  RETRIEVAL  CONCEPT 


When  a  monochromatic  point  object  is  imaged  through  an  optical 
system,  the  observable  is  the  point  spread  function  (PSF)  p(x).  The 
PSF  is  the  modulus  squared  of  the  coherent  system  function  h(x), 

p(x)  =  |h(x)  (  2. 

The  Fourier  transform  of  h(x), 

CO 

H(f)  =  /  h  (x)  exp  (-i2iTfx)  dx , 

— oo 

is  the  coherent  system  function  of  the  optical  system.  Its  modulus, 
A(f) ,  is  a  zero-one  function,  corresponding  to  the  pupil  function, 
and  its  phase  9  (f )  is  proportional  to  the  aberrated  wavefront 
across  the  aperture.  Thus, 

H ( f )  =  A (f 1  exp  lie (f)) . 

Although  we  use  a  one-dimensional  representation  in  this  discussion, 
the  results  are  valid  for  the  physical,  two-dimensional  problem, 
except  as  noted. 

Concisely  stated,  the  phase  retrieval  concept  is  to  estimate 
9(f)  based  on  observation  of  p(x). 

The  basic  algorithm  uses  a  polynomial  expansion  for  9(f), 

00 

9(f)  =  I  ck  9k(f) 
k=l 

where  {0.  (f)}  is  a  suitable  set  of  f-polynomials .  A  phase  estimate 

A 

9(f)  with  a  finite  number  of  c^'s  is  formed  and  an  estimated  OTF , 

P  (f)  ,  is  calculated.  This  is  compared  to  the  measured  OTF,  P(f), 
c 
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Details  on  the  algorithm  that  implements  the  phase  retrieval 
concept  are  given  in  the  final  report  on  the  previous  contract. 

These  details  include  discussions  of  the  search  technique  (a  grid 
search  followed  by  a  modified,  steepest  descent  algorithm  of  the 
Fletcher-Powell  type) ,  the  numerical  procedures  for  calculating  the 
PSF,  the  polynomials  (Zernicke  polynomials,  to  order  24),  and  the 
computer  code  itself.  The  final  report  also  describes  successful 
blind  tests  of  the  algorithm  that  were  conducted  with  Draper  Labs 
and  with  Hughes  Aircraft  Company;  it  describes  a  simple  experimental 
verification  of  the  theory  (point  source,  aberrated  by  a  phase- 
distorted  imaging  system) ;  and  it  gives  some  preliminary  results  on 
extended  objects. 
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ALIASING 


An  array  of  detectors  in  the  focal  plane  spatially  averages 

the  image  over  the  detector  profile  of  each  individual  detector  and 

produces  an  array  of  numbers  that  describes  a  sampled  image.  The 

sampling  interval  is  the  spacing  between  detectors,  measured  from 

center-to-center .  If  the  image  a(x)  is  band-limited  to  a  spatial 

frequency  f  (cycles  per  unit  distance),  then  a(x)  is  completely 

determined  if  spatial  samples  are  taken  at  intervals  smaller  than 

1/(2  f  ),  the  Nyquist  sampling  interval.  In  what  follows  we  will 
max 

set  f  equal  to  1/(XF#),  as  is  the  case  for  a  diffraction-limited 
max 

optical  system. 

For  the  purposes  of  this  study  we  assume  that  the  diode  spacing 
will  be  up  to  five  times  larger  than  the  Nyquist  sampling  interval. 
Thus,  the  image  will  be  undersampled  and  there  will  be  a  severe 
amount  of  "aliasing".  This  aliasing  effect  is  described  as  follows. 
If  the  image  is  a(x)  with  Fourier  transform  A(f) ,  the  sampling 
procedure  produces  a  spectrum  G(f), 

00 

G(f)  =  ^  A(f  _  *fo>  ' 
k=-“> 

where 

xo  =  1//fo 

is  the  spatial  sampling  interval.  If  fQ  is  less  than  twice  the 
highest  frequency  in  A(f),  the  aliases  in  G(f)  overlap  and  the 
spectrum  A(f)  (and,  therefore,  a(x))  cannot  be  recovered.  See 
Figure  6. 


The  effects  of  undersampling  are  more  graphic  when  shown  in 
the  spatial  domain.  Thus,  Figure  7  shows  a  digital  image  of  a  point 
spread  function  sampled  at  the  Nyquist  sampling  rate.  This  PSF  is 
the  result  of  about  one  wave  (peak-to-peak)  of  wavefront  distortion 
in  the  aperture.  In  Figure  8  we  show  another  PSF  with  approximately 
the  same  amount  of  wavefront  distortion.  But  this  has  been  sampled 
by  detectors  that  are  four  times  larger  than  the  Nyquist  sampling 
rate.  Note  the  loss  of  detail. 
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Figure  7.  Part  of  a  32  by  32  point  spread  function 
sampled  at  the  Nyquist  sampling  rate 
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Figure  8.  8  by  8  detected  point  spread  function  undersampled 

by  a  factor  of  4  in  both  X  and  Y  directions 
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This  example  illustrates  the  concern  for  the  utility  of  a  phase 
retrieval  algorithm  when  the  detectors  are  large.  The  image  data 
is  grossly  washed  out.  It  turns  out,  however,  that  the  phase  associ¬ 
ated  with  the  PSF  of  Figure  8  was  accurately  captured  by  the  algorithm, 
even  when  a  small  amount  of  noise  is  added  to  the  observed  PSF. 


In  order  to  study  the  effects  of  aliasing  on  phase  retrieval, 
a  new  computer  program  was  written  which  works  on  PSF's  instead  of 
OTF's.  Square  blocks  of  pixels  may  be  summed  into  one  output  pixel, 
simulating  the  effect  of  using  square  detectors  that  are  an  integral 
number  of  samples  wide.  The  initial,  unaliased  PSF  is  sampled  at 
the  Nyquist  interval  of  \F§/2.  The  pupil  is  a  circular  aperture 
in  a  16  by  16  array,  as  shown  in  Figure  9.  A  phase  aberration  is 
constructed  over  the  aperture  and  the  generalized  pupil  function 
computed.  It  is  buffered  out  to  32  by  32  and  a  Fast  Fourier  Trans¬ 
form  taken,  resulting  in  a  32  by  32  coherent  spread  function.  The 
magnitude  squared  gives  the  PSF. 

4444400000044444 
4440000000000444 
4400000000000044 
4000000000  0  00004 
4000000  0  00  0  00004 
oooooooooooooooo 

OOOOOOOOOOuOOOOO 
OOOOOOOOOOOOOOOO 
OOOOOOOOOOOOOOOO 
OOOOOOOOOOOOOOOO 
OOOOOOOOOOuOOOOO 
4000000000000004 
4000000000  0  00004 
4400000000000044 
4440000000000444 
4444400000044444 


Figure  9.  Pupil  function  used  in  simulations. 

0  indicates  inside  aperture,  4  indicates 
outside  aperture. 
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Simulations  were  run  by  reading  in  an  aberration  (in  terms  of 
Zernicke  polynomials)  and  constructing  a  PSF.  The  detector  model 
was  applied  and  noise  added.  This  is  the  simulated  measured  PSF. 

The  program  attempts  to  deduce  the  aberrations  by  searching  over 
the  parameters  which  are  the  coefficients  of  the  Zernicke  polynomials 
to  find  a  PSF  which  yields  the  best  possible  fit  to  the  measured 
PSF.  The  error  metric  is 

eCp,  =  £  -  Dij)2' 

i.  j 

where  D  is  the  measured  PSF  data  and  P  the  trial  PSF,  which  is  a 
function  of  the  parameters  p. 

The  Zernicke  polynomials  used  are  shown  in  Table  1.  The  results 
for  various  noise  levels  and  amounts  of  aliasing  are  presented  in 
Table  2.  The  row  labelled  "initial"  is  the  aberration  corresponding 
to  the  simulated  measured  PSF.  This  aberration  across  the  aperture 
is  shown  in  Figure  10.  Cases  1  through  7  all  used  detectors  which 
were  averages  of  3  4  by  4  array  of  the  original  Nyquist  sampled 
pixels.  Figure  7  shows  the  Nyquist  sampled  32  by  32  PSF  from  the 
initial  aberration.  Figure  8  shows  the  8  by  8  detected  PSF.  In 
Figure  11  is  the  noisy  PSF  from  case  2. 

Table  1 


Zernicke  Polynomials 

Number 

Form 

Name 

2 

X 

x  tilt 

3 

y 

y  tilt 

4 

2  2 
x  +y 

focus 

5 

2  2 
x  -y 

0°  astigmatism 

6 

2  xy  ^ 

45°  astigmatisi 

7 

x  3) 

x  coma 

8 

y  (  r  2  -  3 ) 

y  coma 

13 


Zernicke  Polynomial  Coefficient 
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Table  2  Results  of  Phase  Retrieval  when  the  PSF  is  undersampled 
by  a  factor  of  4  in  both  X  and  Y  directions 
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Figure  11  Noisy  Detected  Point  Spread  Function 
from  Case  2. 

The  noise  is  added  to  each  pixel  of  the  detected  PSF,  and  is 
independent  and  Gaussian  distributed.  The  standard  deviation  of 
the  noise  is  expressed  in  two  ways  in  Table  2;  as  a  percentage 
of  the  peak  value  of  a  diffraction-limited  PSF  as  detected,  and  of 
the  peak  of  the  aberrated  detected  PSF.  Cases  3  and  4,  and  5  and 
6,  have  the  same  noise  levels  but  have  different  random  realizations 

The  parameter  values  listed  in  each  row  are  the  values  at  which 
the  search  terminated,  and  is  the  final  estimate.  When  the  initial 
parameters  are  corrected  by  the  estimate,  the  residual  phase 
aberration  has  an  RMS  value  across  the  aperture  as  listed  in  the 
column  labelled  "RMS  Phase".  The  Strehl  ratio  of  the  corrected, 
Nvquist  sampled  PSF  is  in  the  column  labelled  "Strehl  Ratio".  Not 
surprisingly,  the  quality  of  the  final  corrected  PSF  deteriorates 
as  the  noise  level  increases,  until  case  7  which  results  in  only 
slioht  improvement. 

From  the  results  in  Table  2  we  see  that  a  PSF,  undersampl.ed  by 
a  factor  of  4,  can  provide  a  good  estimate  of  the  wavefront.  The 
PSF  has  additive  noise  to  a  level  of  about  2%  of  the  diffraction- 
limited  peak. 
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Table  3  shows,  for  the  same  initial  aberration  the  results  of 
increasing  the  detector  size.  Even  with  10  by  10  detectors,  phase 
retrieval  is  still  successful.  Figure  12  shows  the  detected  and  noisy 
detected  PSF's  for  case  13,  with  the  10  by  10  detectors.  Table  3  lists, 
in  the  rows  labelled  "G"  and  "D",  where  the  grid  search  terminated 
and  the  final  result  where  the  Davidon  search  finished. 

Table  4  shows  two  different  aberrations.  Cases  14,  15,  and  16  are 
for  a  detector  size  of  10  by  10.  The  noise-free  case  worked,  and  the 
noisest  case  resulted 'in  significant  improvement,  while  the  less  noisy 
one  essentially  failed  in  the  grid  search  to  find  a  good  starting  point 
for  the  Davidon  search.  Cases  17  and  18  show  another  aberration  failing 
with  5  by  5  detectors,  both  for  noise-free  and  noisy  data.  With  2  by  2 
or  4  by  4  detectors,  phase  retrieval  was  successful. 

These  results  show  that  aliasing  is  not  a  major  problem  to  phase 
retrieval.  There  are  still  problems  with  regard  to  doing  a  fine 
enough  mesh  on  the  grid  search  to  insure  detection  of  the  global 
minimum,  but  aliasing  does  not  seem  to  create  any  new  difficulties. 


Figure  12 
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Zernicke  Polynomial  Coef 


Results  of  Phase  Retrieval 


4  NOTES  ON  EFFICIENT  SEARCH  TECHNIQUES 

Previously  proposed  phase  retrieval  search  methods  have  involved 
comparing  a  measured  PSF  (or,  equivalently,  an  OTF)  with  a  large 
number  of  PSF’s  from  a  grid  of  parameter  values  covering  the  region 
of  parameter  space  which  is  expected  to  include  the  actual  aberration. 
These  PSF '  s  have  always  been  recomputed  for  each  search,  even  though 
they  do  not  change.  A  faster  method,  which  actually  uses  less 
computational  hardware,  is  to  compute  the  PSF's  once  and  store  them 
on  disk,  and  retrieve  them  as  needed  during  the  search. 

To  work  out  the  speed  of  this  operation,  suppose  that  we  wish 
to  examine,  in  the  grid  search,  three  values  for  each  of  thirteen 
parameters.  This  yields  1.6  x  10^  PSF's.  (Five  values  for  nine 
parameters  gives  a  slightly  larger  number.)  Storing  an  eight  by 
eight  PSF,  quantized  to  256  levels  (so  that  each  pixel  will  fit  ir 
a  single  byte)  for  each  grid  point  would  take  102  megabytes,  which 
is  within  the  capacity  of  a  single  large  disk  drive.  If  the  data  is 
read  off  the  disk  in  an  optimum  manner  so  as  not  to  waste  disk  head 
movements,  three  minutes  are  taken  to  transfer  all  the  data  to  the 
CPU.  To  compute  the  error  metric  on  a  64  pixel  PSF  takes  128 

g 

additions  and  64  multiplications.  The  total  for  all  1.5  x  10  PSF's 
6  6 

is  192  •  10  additions  and  96  *  10  multiplications .  To  perform  the 
computation  in  three  minutes  of  disk  read  time  requires  a  speed  of 
500  nsec  per  addition  and  1000  nsec  per  multiplication.  This  is 
within  the  capability  of  a  PDP  11/70  doing  integer  arithmetic. 
(Floating  point  is  not  needed  for  these  calculations.)  Since  the 
CPU  would  be  active  at  the  same  time  as  the  disk  transfer  were  takinq 
place,  the  computation  time  and  I/O  time  could  be  totally  overlapped 
so  that  only  three  minutes  are  required  for  the  entire  search. 
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There  remain.--  the  problem  of  the  second  grid  search.  This 
is  performed  at  half  the  grid  spacing  of  the  first  search  in  the 
region  immediately  around  the  best  points  of  the  first  search. 

It  has  proven  necessary  only  occasionally  to  actually  do  the  second 
grid  search  in  the  simulations  that  have  been  run  previously,  but 
for  higher  dimensions  of  parameter  space  the  distance  from  a  point 
in  the  center  of  any  hypercube  formed  by  the  grid  points  grows  as 
the  square  root  of  the  dimensionality.  Prestcring  PSF's  calculated 
for  parameter  values  on  a  finer  grid  mesh  would  take  far  too  much 
storage.  Calculating  them  as  needed  for  a  second  grid  search  would 
dwarf  the  time  of  the  first  search.  If  the  first  grid  spacing  is 
small  enough,  then  the  second  search  can  be  avoided.  Alternatively, 
a  downhill  search  can  be  started  from  each  of  the  best  survivors  of 
the  first  search.  The  downhill  search  time  rises  only  linearly 
with  the  dimensionality. 

A  faster  method  for  the  first  grid  search  involves  the  use  of 
moments,  or  some  other  characterization  of  distinctive  properties 
of  the  PSF.  Five  moment?,  as  we  use  in  Section  5,  seem  a  good  number 
The  PSF's  used  for  the  first  grid  search  are  sorted  on  their  moments 
into  15  bins  adjusted  for  each  moment  to  span  the  range  of  moment 
values,  and  perhaps  so  that  approximately  equal  numbers  of  PSF's 
fall  into  each  bin.  The  PSF's  are  placed  on  the  disk  in  an  order 
beginning  with  the  PSF's  where  moments  all  fall  into  the  first  bins, 
stepping  through  all  bins  on  the  fifth  moment,  then  moving  onto  the 
second  bin  of  the  fourth  moment  and  stepping  through  all  fifth 
moment  bins,  and  cortinuing  until  the  PSF's  in  the  fifteenth  bin 
for  all  moments  are  written.  This  order  is  listed  in  Table  5. 


TABLE  5 

Order  of  PSF's  Sorted  in  Moments 


1 

1 


1 

1 


M3  M4 

1  1 

1  1 


1 

2 


1  1  1  1  15 

1112  1 


1  1  1  2  15 


1  1  1  15  15 

112  11 

15  15  15  15  15 


Thxs  is  a  total  of  15"  ,  or  759,375  possible  slots  for  a  PSF  to  be 
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assigned  to.  If  we  have  3  grid  points,  or  1,594,323,  there  will 
be  an  average  of  two  PSF's  in  each  slot,  although  considerable  varia¬ 
tion  may  be  expected.  Given  a  measured  PSF,  its  five  moments  may 
be  calculated  and  a  decision  made  as  to  which  bin  each  moment  falls 
in,  in  negligible  time.  A  way  of  pointing  to  the  area  on  disk  which 
holds  the  PSF's  corresponding  to  the  five  moments  is  needed.  Such 
a  table,  with  759,375  entries  may  be  constructed.  Its  size  is  such 
that  it  too  must  be  on  disk.  Each  table  entry  holds  an  address  on 
the  disk  for  where  PSF  data  begins  for  that  particular  slot.  Data 
continues  on  the  disk  until  the  next  address  in  the  table  is  reached. 
Access  into  table  is  easily  computable;  if  the  bin  values  for  each 

moment  are  i,  j,  k,  C,  and  m,  then  element  number  (i— 1) *15^  +  (j-l)*15 
2 

+  (k— 1)  *15  +  (?.-l)*15  +  m  is  the  proper  table  address.  Generating 
the  sorted  PSF  file  and  the  address  table  will  require  a  great  deal 
of  computer  time,  but  it  only  needs  to  be  done  once. 
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Examining  only  the  bins  into  which  the  moments  of  measured  PSF 
fall  is  unwise,  since  noise  and  the  presence  of  other  aberrations 
will  perturb  the  moments.  A  safer  course  is  to  examine  *-hree  bins 
for  each  moment,  the  one  into  which  the  measured  moment  fails  and 
the  ones  on  either  side.  This  corresponds  to  20%  of  the  range  of 
moment  values  for  each  moment.  For  five  moments,  only  .00032  of 
the  full  set  of  grid  points  survive.  Taking  three  bins  for  each 
moment  results  in  making  3^,  or  243,  references  in  the  address  table, 
which  will  point  to,  on  the  average,  a  total  of  486  PSF's.  Each  of 
these  PSF's  may  be  compared  with  the  measured  PSF  and  a  metric 
computed.  This  should  take  less  than  30  seconds  to  do.  Downhill 
searches  can  be  started  from  the  parameter  sets  associated  with  the 
PSF's  with  the  lowest  metrics.  The  parameters  from  the  search  that 
terminates  at  the  smallest  metric  are  the  final  estimate. 
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THE  USE  OF  MOMENTS  TO  SPEED  UP  THE  ALGORITHM 


The  search  method  which  has  been  used  for  phase  retrieval 
applications  consists  of  an  initial  search  over  a  rectangular  grid 
of  points  in  parameter  space,  a  second  search  around  the  best  (lowest 
metric)  grid  points  at  half  the  first  grid  spacing,  and  a  downhill 
Davidon  search  from  the  best  point  of  the  second  grid  search.  The 
initial  grid  search  is  the  most  crucial  to  obtaining  a  good  phase 
estimate,  and  is  quite  time-consuming  to  compute  the  metric  of  the 
sum  of  the  squares  of  the  differences  between  the  measured  and  trial 
PSF's  point  by  point. 

A  method  was  sought  that  would  allow,  by  means  of  a  quick 
comparison  of  only  a  very  few  numbers  associated  with  each  PSF, 
whether  or  not  a  trial  PSF  is  worthy  of  further  consideration.  '  A 
natural  approach  is  to  use  moments  of  the  PSF  distribution.  The 
most  fundamental  moment  is  the  center  of  gravity,  defined  by 

E  =  J P{x,y)dxdy 

x  =  f x  P(x,y)dxdy/E 

c  .  cj  •  * 

yc.g.  =  /y  p(x/y)dxdy/E- 

Given  that  the  ideal  position  of  the  point  object  is  known,  the 
location  of  the  c.g.  (in  the  absence  of  noise)  exactly  determines 
the  tilt  parameters  of  the  aberration.  A  slightly  modified  set  of 
Zernicke  polynomials  may  be  constructed  such  that  only  the  tilt 
terms  move  the  c.g. 
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Previous  work  has  shown  the  c.g.  to  be  rather  sensitive  to 
quantization  noise  (and  presumably  any  other  sort  of  noise)  on  the 
PSF.  In  the  work  reported  here,  we  have  assumed  that  the  tilt 
problem  will  be  dealt  with  in  some  appropriate  fashion.  The 
aberrations  tested  had  no  tilt  errors. 

Five  more  moments  around  the  c.g.  were  defined;  the  weighting 
factors  were  taken  to  be  the  aberration  polynomials  themselves.  They 
are: 


s 

I—1 

II 

J(x’ 2+y*  2)  P  dxdy/E 

m2  = 

J  x'  y'  P  dxdy/E 

M3  = 

/ (x,2-y'2)  P  dxdy/E 

M4  = 

J*  x'r  P  dxdy/E 

MS  = 

J  y' r2  P  dxdy/E  , 

!  X' 

=  X  -  X 

eg 

y 

ii 

i 

in 

Some  intuitive  meanings  may  be  assigned  to  these  moments.  is  an 

indication  of  the  overall  width  of  the  PSF.  M2  is  large  and  positive 
for  a  PSF  of  an  elliptical  sort  of  shape,  with  the  long  axis  oriented 
at  45°  to  the  X  axis.  It  is  negative  for  the  long  axis  along  135°. 

The  greater  the  asymmetry  between  the  long  and  short  axes,  the  larger 
the  magnitude  of  M2-  is  similar  to  M2  except  that  the  orientations 

for  maximum  magnitude  are  0°  and  90°.  Orientations  of  the  long  axis 
at  other  angles  give  contributions  to  both  M2  and  and 

measure  the  asymmetry  along  either  the  X  or  Y  directions  relative 
to  the  c.g.  A  PSF  like  Figure  13  would  have  a  positive  value. 
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Center  of  Gravity 


Figure  13  PSF  with  large  positive 

Sofware  was  written  to  calculate  these  moments  from  sim'  lated 
PSF's.  An  investigation  was  made  of  the  sensitivity  of  the  moments 
to  noise  on  the  PSF.  Table  6a  shows  the  five  moments  and  the  c.g. 
for  the  PSF  from  an  aberration  consisting  of  1/4  wave  of  45° 
astigmatism  and  1/4  wave  of  x-coma.  The  sample  spacing  of  the  PSF 
is  \F#/2.  The  first  line  shows  results  with  no  no  ;e.  The  next  two 
lines  show  the  averages  and  standard  deviations  of  the  moments  for 
100  realizations  of  signal-independent  white  noise  on  the  PSF.  The 
noise  variance  was  1.85s  of  the  maximum  value  of  the  diffraction- 
limited  PSF.  Since  this  aberration  has  a  Strehl  ratio  of  .468,  the 
variance  was  3. OS’s  of  the  peak  value  of  the  aberrated  PSF.  The 
variation  of  the  moments  is  too  large  for  them  to  be  useful.  Figures 
14  and  15  show  this  PSF  noise-free  and  with  noise. 
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Moments  and  Centers  of  Gravity 
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To  reduce  the  effects  of  noise,  a  Gaussian  window  was  applied 
to  the  moment  calculations,  as  in 

E  =  /e“(x'  +y'  )  /w  P  dxdy/E 

Mx  =  /  (x'2+y'2)  e_(x  ,+y  1 /w  P  dxdy/E 

The  width  w  was  taken  to  be  5  sample  spaces.  This  eliminates  the 
contribution  from  small  noise  fluctuations  located  at  large  distances 
from  the  c.g.,  which  weight  heavily  into  the  moments.  This  significantly 
reduced  the  variations,  but  they  are  still  too  large,  as  shown  in 
Table  6b. 

The  successful  strategy  was  to  also  apply  the  same  Gaussian 
window  to  the  c.g.  calculation.  The  c.g.  was  first  computed  as 
before,  and  then  recomputed  with  the  Gaussian  window  centered  on 
the  first  c.g.  estimate.  As  may  be  seen  in  Table  6c.  the  standard 
deviation  of  the  c.g.,  and  consequently  the  five  moments,  were 
reduced  by  factors  of  four. 

The  following  procedure  was  used  to  test  the  sensitivity  of 
a  search  on  the  moments  to  noise  and  aliasing.  A  table  of  moments 
was  generated  for  1575  PSF's  on  the  grid  of  five  values  of  each 
aberration  coefficient  (-.  5  ,  -.  25  ,  0  .,  .  25  ,  and  .5  waves)  for  five 
different  coefficients  (defocus,  0°  astigmatism,  45°  astigmatism, 
x  coma,  and  y  coma).  (The  number  of  such  grid  points  is  reduced 
from  55  by  a  factor  of  nearly  one  half  due  to  symmetry.)  A 
particular  aberration  was  chosen  (.25  defocus,  .5  0°  astigmatism, 

-.25  y-coma) ,  which  was  exactly  one  of  the  grid  points.  Moments 
were  calculated  from  the  corresponding  PSF's  for  various  amounts 
of  noise  and  aliasing.  The  table  of  moments  was  searched  for  the 
best  fit  to  the  moments  from  the  noisy  PSF's.  In  the  absence  of 
noise,  the  moments  in  the  table  corresponding  to  the  same  aberration 
as  the  test  case  would  be  a  perfect  match.  With  enough  noise,  some 
other  aberrations  will  have  better  fits  to  the  trial  moments. 
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At  each  noise  level,  ona  hundred  different  noise  realizations 
on  the  PSF  were  run  through  the  moment  calculation  and  search. 

Counts  were  made  of  the  number  of  cases  in  which  more  than  a  certain 
number  of  aberrations  had  better  fits  to  the  noisy  trial  moments 
than  the  mo...^nts  corresponding  to  the  actual  aberration.  The  worst 
case  (with  the  largest  number  better  than  the  actual  aberration) 
was  also  noted.  The  v^orst  case  indicates  how  deep  in  the  list  of 
best  fits  to  the  moments  a  search  over  the  PSF  itself  must  go, 
to  insure  inclusion  of  the  correct  PSF  in  the  PSF's  examined. 

Table  7  presents  the  results,  in  terms  of  the  rank  of  the  correct 
grid  point  (correct  phase) . 

From  Table  7  we  see  that,  for  example,  with  a  detector  size 
equal  to  the  Nyquist  sampling  interval  (lxl  case)  and  with  50% 
dependent  noise,  there  were  13  noisy  PSF's  (out  of  100  tested)  where 
the  correct  phase  was  not  ranked  number  1.  The  worst  rank  was 
number  10. 

The  noise  was  added  to  each  point  of  the  PSF  by 

PSF'(X,Y)  =  PSF  (X ,  Y)  (1  +  N1-A[))  + 

where  and  ^  are  independent  random  variables  with  zero  mean 
and  a  variance  of  one.  AQ  and  A^.  are  the  amplitudes  of  the  dependent 
and  independent  noise  terms.  The  independent  noise  amplitude  is 
given  as  an  absolute  number  and  as  a  percentage  of  the  peak  value 
of  a  diffraction-limited  PSF  and  of  the  peak  value  of  the  aberrated 
PSF. 

In  the  case  of  no  aliasing,  a  remarkable  amount  of  signal 
dependent  noise  can  be  tolerated.  Even  with  A^  equal  to  70%,  the 
worst  case  had  only  28  grid  points  with  a  better  fit  to  the  moments. 
The  permissible  dependent  noise  drops  roughly  with  the  linear  amount 
of  aliasing.  The  number  of  PSF  data  points  drops  with  the  square 
of  the  aliasing  factor.  The  permissible  signal  independent  noise 
expressed  as  a  percentage  of  the  peak  value  of  the  aberrated  PSF 
also  drops  in  a  linear  manner,  but  a  much  smaller  amount  of  noise 
can  be  tolerated. 
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Table  7 


Cumulative  Histogram  of  the  Rank  of 
the  Correct  Grid  Point 

For  each  noise  level  100  noisy  PSF’s  are  subjected  to  the  grid 
search. 

K  =  Index  of  the  Kth  noisy  PSF,  K  =  1,2,  .  .  100, 

X  =  Rank  of  the  correct  grid  point  out  of  all  1575  points 

*»  t  h 

tested  in  the  K-  search. 


6  PHASE  DIVERSITY  IMAGING 

In  the  previous  contract  we  developed  several  algorithms  to 
perform  phase  retrieval  on  distorted  images  of  extended  objects. 

The  optical  system  introduces  a  phase  aberration  0(f)  with  its 
corresponding  PSF  p(x).  The  observed  image  is 

z(x)  =  o(x)  *  p  ( x )  +  n(x), 

where  means  convolution,  o (x)  is  the  unknown,  extended  object 

and  n(x)  is  noise. 

Since  o(x)  is  unknown,  we  cannot  use  the  algorithm  of  Section  2 
to  find  8(f).  However,  phase  retrieval  is  still  possible  if  we  can 
acquire  two  images  z^  (x)  and  z,,(x)  of  the  object.  These  two  images 
allow  us  to  make  a  joint  estimate  of  o(x)  and  0(f).  The  details  are 
in  the  previous  Final  Report  (Appendix  C) .  Briefly,  for  an  estimated 

A  A  A 

object  o(x)  and  estimated  PSF's  p-^(x)  and  (x) ,  we  calculate  an 
error  metric 

E  =  f[  z1(x)  -  ofx^p^x)]2  dx 
f{ z2(x)  -  o(x)*p2(x)]2  dx. 

When  we  use  Parseval's  theorem  to  perform  this  calculation  in  the 
spatial  frequency  domain,  we  can  minimize  E  by  choice  of  0(f):  namely 

P*  (f) Z,  (f)  +  P* (f) Z9  (f) 

0(f)  =  — - , - 

I  px (f ) I  +  !p2(£) I 

We  substitute  this  0(f)  into  the  E-metric  calculation  and  find  the 
0(f)  that  further  minimizes  the  integrals.  The  0-search  is  now 
identical  to  that  of  Section  2. 

To  implement  this  algorithm  we  have  chosen  to  add  a  known 
phase  0(f)  to  the  already  phase-aberrated  optical  system.  With 
0(f)  =  0  we  obtain  the  first  image  z^(x);  with  p(f)  ?  0  we 
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Now  we  show  a  computer  simulation  of  the  technique.  A  point 
is  imaged  through  an  aberrated  optical  system,  with  and  without  an 
intentional  defocus.  The  two  point  spread  functions  are  measured 
and  are  used  to  estimate  the  object.  The  estimation  procedure,  of 
course,  has  no  knowledge  about  the  object. 

In  Figure  17  we  show  a  two-dimensional  phase  aberration  8  in 
a  circular,  unapodized  aperture.  The  resulting  point  spread  func¬ 
tion  of  the  optical  system  is  shown  in  Figure  18.  Next,  we  inten¬ 
tionally  defocus  the  system  (one-half  a  wave  of  quadratic  phase  from 
the  center  to  the  edge  of  the  aperture)  .  The  new  phase  (0  +  <j>)  is 
shown  in  Figure  19  and  the  corresponding  point  spread  function  in  20. 

In  this  example  we  add  no  noise  so  that  z-^  and  z2,  the  observed 
signals,  are  those  shown  in  Figures  18  and  20.  Their  Fourier  trans¬ 
forms,  and  Z2,  are  the  inputs  to  the  block  diagram  of  Figure  16. 
The  outputs  are  0  (an  estimate  of  0)  and  6  (an  estimate  of  the 
object's  spectrum).  The  former  is  shown  in  Figure  21  and  should  be 
compared  with  Figure  18.  The  agreement  is  excellent.  The  inverse 
Fourier  transform  of  6  is  the  estimated  object.  This  is  shown  in 
Figure  22. 


35 


Ol.fltt 


0 . 000 


SSS  S8SSS&SSES 

e  in  ©  m  s  in  s  s  in  ©  m  s 
—  — *  r\.  r\j  m  on  •■=*  m  m  ld  —■  r'-- 


G«i  +  XO«-KN>-XXM-ft  i  _ 


f 


mu 

[MU 

t\\\\  A  v  \  \  \ 

iWAVa '  \  \  v 


v\A'V'  ’ 

ft V\ 'AH 


WXW' 


'V  A 


\v\  j- 

\\  \  \  ! 


<s^\\ 

Nvv\\ 
\\  \  \  > 

Any- 

av ) 


/ 1 , 

(J\ 

y  ^ 


Next  we  present  a  computer  simulation  where  the  object  is  a 
segment  of  a  water  tower  as  shown  in  Figure  23(a).  The  initial 
system  aberration  is  the  same  S  as  in  Figure  17.  The  resulting 
image  is  shown  in  Figure  23(b)  and  the  defocused  image  is  shown  in 
Figure  23(c).  The  Fourier  transforms  of  the  images  in  23(b)  and 
23(c)  are  the  inputs  to  the  algorithms  of  Figure  16.  The  result  is 
shown  in  Figure  23(d). 


Figure  23. 
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When  the  phase  diversity  is  quadratic,  successive  images 
represent  the  images  taken  at  different  values  of  defocus.  This 
is  common  practice  in  electron  microscopy  for  reconstruction  of  the 
phase  of  an  object  under  observation.  It  seems  reasonable  to  expect 
that  a  particular  kind  of  phase  diversity,  other  than  quadratic, 
will  allow  better  object  reconstruction  for  particular  kinds  of 
channel  distortions. 

The  main  utility  of  the  technique  may  be  the  fact  that  the 
adaptation  is  not  done  in  closed  loop  fashion.  Thus,  the  adaptive 
optics  does  not  have  to  track  the  temporal  variation  in  the  channel. 
The  phase  diversity  allows  the  images  to  be  unscrambled  in  post¬ 
processing.  On  the  other  hand,  it  would  probably  be  advantageous 
to. include  a  closed  loop  with  a  long  time  constant  to  remove  gross 
distortions  in  the  channel.  This  concept  could  greatly  reduce  the 
burden  on  the  adaptive  control  mechanism  while  increasing  the  burden 
on  post-processing  hardware. 
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LABORATORY  EXPERIMENT 


7 . 1  Laboratory  Set-Up 

A  procedure  has  been  developed  to  experimentally  test  the  phase 
retrieval  concept  for  extended  objects.  This  procedure  is  similar 
to  earlier  experiments  in  which  image  data  from  an  optical  system 
is  recorded.  The  optical  system  contains  a  known  phase  aberration 
function.  This  function  is  contained  on  a  clear  substrate  and  is 
separable  from  the  rest  of  the  optical  system.  Analysis  of  the 
aberration  alone  is  possible  by  interferometric  techniques.  This 
result  provides  the  knowledge  of  the  function  against  which  the 
retrieved  phase  is  to  be  compared.  To  supply  sufficient  data  for 
phase  retrieval  of  the  aberration  function  in  this  experiment,  the 
image  is  recorded  at  two  positions  of  known  focus  difference. 

The  imaging  system  is  sketched  in  Figure  24.  It  is  a  1:1 
projector  with  condenser  illumination.  The  lamp  illuminates  a  ground 
glass  diffuser.  The  diffused  lamp  image  is  relayed  by  the  condenser 
to  the  system  aperture  stop.  At  the  condenser  are  filters  to  limit 
the  spectral  transmittance  of  the  system  to  red  light.  The  aperture 
stop  is  located  midway  between  two  acromats.  Each  acromat  is  260mm 
in  focal  length.  A  1:1  objective  is  made  by  placing  the  object  and 
image  conjugate  at  the  focus  of  each  lens.  The  region  between  lenses 
contains  collimated  rays.  It  is  in  this  region  that  the  phase 
aberration  is  inserted.  The  effect  of  the  aberration  on  the  optical 
system  is  therefore  the  same  as  its  effect  in  the  collimated  beam 
of  an  interferometer.  This  maximizes  the  agreement  between  the  inter¬ 
ferometric  measure  of  the  aberration  and  that  found  by  phase  retrieval. 

doublet 
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The  second  doublet  focuses  the  image  onto  a  plane  which  is 
scanned  by  a  photodiode  array.  The  array  has  1024  elements  with 
15pm  center  to  center  spacing.  The  array  is  stopped  15pm  at  a  time 
to  produce  a  sampled  image  of  1024x1024  pixels. 

The  phase  aberration  can  be  inserted  or  removed  repeatably 
between  the  doublets.  A  singlet  of  1000mm  focal  length  can  also 
be  inserted  behind  the  second  lens.  This  lens  controllably  defocuses 
the  image . 

The  spectral  distribution  of  the  illumination  was  chosen  to  be 
primarily  red,  the  region  of  peak  diode  sensitivity.  It  was  desired 
to  be  a  broad  bandwidth  for  incoherent  imaging  yet  narrow  enough 
to  use  its  central  wavelength  for  calculations  which  are  wavelength 
dependent.  The  relative  spectral  response  of  the  system  is  centered 
at  640nm  with  a  spread,  roughly,  of  lOOnm.  It  is  the  combination  of 
the  lamp  distribution  at  2850K,  an  infrared  sharp  cutoff  filter,  a 
red  broadband  interference  filter  and  the  spectral  response  of  the 
photodiodes . 

The  aperture  size  of  2.2mm  and  focal  distance  of  260mm  gives 
an  F#  of  118.  This  was  selected  to  give  a  diffraction  limited 
spread  function  width  of  10  diodes  (150pm) ,  a  suitable  sampling 
density  for  the  phase  retrieval  algorithm. 

A  two  beam  interferometer  was  constructed  to  provide  a  direct 
measure  of  the  aberrations  used.  The  setup  was  standard.  One  beam 
of  parallel  wavefronts  is  the  reference  and  the  other  is  the  signal 
beam.  The  aberration  is  placed  in  the  signal  beam  and  the  two 
beams,  when  recombined,  display  an  interference  pattern.  This 
pattern  is  recorded  on  film  as  an  interf erogram.  An  example  is 
shown  in  Figure  25.  This  aberration  appears  to  have  about  two  waves 
of  distortion  and  was  used  to  produce  the  images  shown  later. 
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7.2  Data 


Figures  26-29  show  a  series  of  images.  These  represent  the 
four  imaging  conditions  which  were  available  in  this  experiment. 
Figure  26  shows  an  image  of  a  section  of  a  tri-bar  target  produced 
by  the  diffraction  limited  system.  Figure  27  shows  the  same 
image  with  +0.9  waves  of  defocus  added  by  including  the  defocusing 
lens  in  the  system.  Figure  28  shows  the  system  with  the  aberration 
of  Figure  25  added  at  the  aperture  stop.  Figure  29  shows  the  effect 
of  the  same  aberration  and  +0.9  waves  of  defocus. 

In  later  data  processing  we  used  images  in  Figures  28  and  29 
to  determine,  by  our  two-focal-plane  algorithm,  the  wavefront 
distortion . 
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7.3  Determination  of  the  Quadratic  Phase 


Two  measured  diffraction  limited  images  (one  with  quadratic 
phase  added  to  the  phase  and  one  without)  were  input  to  the  phase- 
retrieval  algorithm.  The  images  used  are  given  in  the  preceding 
section  (Figures  26  and  27).  The  idea  was  to  determine,  from  these 
two  images,  the  estimated  phase  distortion  across  the  aperture  of 
the  system.  By  imposing  the  known  focus  error  between  the  two 
images,  the  calculated  phase-estimate  should  be  a  constant  and  thus 
verify  that,  the  estimated  focus  error  is  correct.  The  first  step 
in  processing  the  data  was  to  obtain  a  smaller  region  of  the  images 
to  use  as  input  to  the  algorithm,  for  the  simple  reason  that  the 
entire  image  would  have  been  burdensome  computationally.  Corresponding 
80x80  pixel  regions  of  each  image  were  chosen  for  processing.  The 
actual  region  corresponds  to  the  area  surrounding  the  number  "4"  in 
the  right  side  of  the  tri-bar  target  in  Figure  26. 

The  first  moment  (center  of  gravity)  of  each  image  was  calculated 
so  that  each  image  could  be  aligned  to  remove  major  tilts  in  the  phase 
estimate.  The  slow  Fourier  transform  was  performed  on  each  section, 
each  output  transform  being  a  17x17  sampled  OTF.  The  transform  was 
sampled  so  that  the  outer  point  corresponds  to  the  diffraction-limited 
cutoff  of  the  system.  These  two  OTF's  were  then  input  to  the  two- 
focal  plane  phase-retrieval  search  algorithm  in  order  to  estimate  the 
phase.  The  input  focus  error  was  assumed  to  be  0.9  waves. 

The  phase  estimate  found  from  the  search  algorithm  is  given  in 
Figure  30.  The  RMS  phase  was  calculated  to  be  .022  waves  and  is  .06 
waves  peak  to  peak.  The  result  is  consistent  with  a  diffraction- 
limited  phase  to  within  the  errors  expected  in  the  experimental 
procedure. 


7 . 4  Determination  of  the  Unknown  Aberration 

The  two  aberrated  images  {one  with  focus  error  added  and  one 
without)  were  input  to  the  phase  retrieval  algorithm.  The  same 
corresponding  80x80  images  section  used  for  the  diffraction  limited 
(Number  "4")  was  processed  using  the  exact  same  procedure  as  the 
diffraction  limited  case.  The  phase  retrieval  algorithm  was  allowed 
to  search  over  20  Zernike  coefficients  in  order  to  fit  the  correct 
phase.  The  Zernike  coefficients  are  given  in  Figure  31  and  a  contour 
plot  of  the  phase 
that  was  measured 


is  given  in  Figure  32.  The  actual  phase  aberration 
interferometrically  is  given  in  Figure  33. 
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phase  retrieval  algorithm  -  applied  to  two-focus 
problem  (aberrated  case) 
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BLIND  TEST 


To  test  the  phase  retrieval  algorithm,  the  customer  requested 
that  personnel  from  Draper  Lab  give  us  several  noisy,  undersampled 
point  spread  functions.  These  were  prepared  by  Dr.  V.N.  Mahajan 
and  his  staff  at  the  Optical  Systems  Division  in  Cambridge,  MA. 

We  subsequently  received  three  aberrated  point  spread  functions, 
each  "sampled  by  an  8x8  array  of  detector  elements.  The  width  of 
the  detector  elements  is  2.13XF,  where  X  is  the  wavelength  of  the 
object  radiation  and  F  is  the  focal  ratio  (f #)  of  the  optical 
system."*  The  three  point  spread  functions  are  shown  in  Figure  34. 
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Figure  34.  PSF’s  for  the  blind  test 


* 

Letter  of  June  30,  1981  to  R.  Gonsalves  from  V.  Mahajan 
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We  exerc ' sed  the  phase  retrieval  algorithm  on  the  three  point 
spread  functions  (PSF)  as  follows.  First,  we  left  the  sampling 
interval  at  2.13AF#,  even  though  our  algorithm  assumes  the  prefix 
constant  to  be  a  integer  (we  used  2) .  For  each  PSF  we  assumed  nat 
the  unknown  phase  could  be  approximated  by  8,  by  15,  and  by  23 
terms  in  a  Zernike  polynomial  expansion.  The  8  and  15  term  expansions 
were  determined  with  the  all-zero  polynomial  as  a  starting  point  in 
our  search.  The  23  term  expansion  used  the  8  term  solution  as  a 
starting  point. 

Next  we  interpolated  the  PSF's  to  give  us  an  8x8  array  sampled 
at  2AF# .  The  searches  described  above  were  repeated. 

In  Figure  35  we  list  the  Zernike  polynomials. 

Figure  36  is  a  summary  of  the  data  processing  for  the  18  cases. 

Figures  37  through  82  give  the  results.  For  example,  Figure  37 
shows  the  input  PSF,  the  estimated  PSF,  the  mean  square  error  between 
them,  and  the  coefficients  of  the  Zernike  polynomials  for  the 
estimated  phase.  Figure  38  shows  the  estimated  phase  itself.  This 
first  pair  of  figures  is  for  case  1,  sampling  interval  =  2.13XF#, 
and  23  Zernike  polynomials. 
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Figure  36.  Data  Processing  Summary 
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Figure  42.  Blind  test  results  for  case  1, 
2.13XF#,  8  polynomials 
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Figure  43.  Blind  test  results  for  case  2, 
2.13AF#,  23  polynomials 
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Figure  44.  Blind  test  results  for  case  2, 
2.13XF#,  23  polynomials 
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Figure  45.  Blind  test  results  for  case  2, 
2.13XF#,  15  polynomials 
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Figure  46.  Blind  test  results  for  case  2, 
2.13AF#,  15  polynomials 
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Figure  47.  Blind  test  results  for  case  2, 
2.13XF#,  8  polynomials 
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Figure  48.  Blind  test  results  for  case  2, 
2.13AF#,  8  polynomials 


73 


I 


PHRMN25.LOG;  1  28-JUL-198t  15:38:26.96 

NUMBER  OF  ZERNIKE  POLYNOMIALS  =  23 


NUMBER  OF 
ACTUAL  NUMBER  OF 

INPUT  PSF 


34 

-43 

14 

87 

- 1  6 

-21 

51 

108 

30 

41 

22 

-9 

27 

22 

43 

-7 

-5 

53 

-81 

-32 

166 

-27 

24 

-25 

53 

-63 

7 

13 

375 

29 

13 

40 

57 

-2 

71 

15 

122 

-2 

14 

69 

16 

21 

-6 

-60 

56 

-41 

16 

-14 

8 

-37 

12 

-2 

-14 

-94 

-91 

-40 

34 

-17 

-58 

45 

-49 

-41 

-23 

-70 

MEAN  SQUARE  ERROR 


COEFFICIENTS  OF  THE 
ZERNIKE  POLYNOMIALS 


ITERATIONS  =100 
ITERATIONS  =  21 


ESTIMATED  PSF 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

2 

5 

1 

0 

0 

0 

0 

0 

1 

139 

2 

1 

0 

0 

0 

1 

14 

348 

20 

3 

0 

0 

0 

1 

4 

95 

2 

1 

1 

0 

0 

0 

1 

1 1 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

*  113471.9531 


K 

PCK) 

t 

0.00000 

2 

-0. 72504 

3 

0.33835 

4 

0.14078 

5 

0.30064 

6 

-0.01004 

7 

0.04718 

8 

0.02084 

9 

0.02034 

10 

0.01474 

11 

-0.00880 

12 

0.00405 

13 

0.01477 

14 

0.01639 

15 

0.03228 

16 

-0.00119 

17 

0.00482 

18 

-0.01799 

19 

-0.03260 

20 

0.04583 

21 

-0.01725 

22 

0.00017 

23 

-0.00588 

Figure  49.  Blind  test  results  for  case  3, 
2.13AF#,  23  polynomials 
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Figure  50.  Blind  test  results  for  case  3, 
2.13XF#,  23  polynomials 
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Figure  51.  Blind  test  results  for  case  3, 
2.13AF#,  15  polynomials 
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Figure  62.  Blind  test  results  for  case  3, 
2.13XF#,  15  polynomials 
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Figure  53.  Blind  test  results  for  case  3, 
2.13XF#,  8  polynomials 
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Figure  54.  Blind  test  results  for  case  3, 
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Blind  test  results  for  case  1, 
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Figure  58.  Blind  test  results  for  case  1, 
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Figure  60.  Blind  test  results  for  case  1 
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Figure  62.  Blind  test  results  for  case  2, 
2AF #,  23  polynomials 
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Figure  63.  Blind  test  results  for  case  2, 
2XF# ,  15  polynomials 
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Figure  64.  Blind  test  results  for  case  2, 
2 AF# ,  15  polynomials 
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Figure  65.  Blind  test  results  for  case  2, 
2  AF# ,  8  polynomials 
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Figure  66, 


Blind  test  results  for  case  2, 
2XF# ,  8  polynomials 
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Figure  67.  Blind  test  results  for  case  3, 
2AF#,  23  polynomials 
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Figure  68.  Blind  test  results  for  case  3, 
2XF#,  23  polynomials 
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Figure  69.  Blind  test  results  for  case  3, 
2AF#f  15  polynomials 
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Figure  70.  Blind  test  results  for  case  3, 
2XF#,  15  polynomials 
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Figure  71.  Blind  test  results  for  case  1, 
2XF#,  8  polynomials 
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Figure  72.  Blind  test  results  for  case  3, 
2AF#,  8  polynomials 


Notes  Added  in  Revision 
January  27,  1982 


The  results  of  our  phase  estimates,  pages  62  to  97,  were  given 
to  V.  Mahajan  of  Draper  Labs  on  August  14,  1981.  Several  problems 
with  our  data  were  identified  by  Draper  personnel  when  they  tried  to 
reconcile  our  phases  with  their  input  phases.  Our  Zernike  polynomials 
are  unnormalized  (the  search  does  not  require  normalized  polynomials) , 
there  is  an  apparent  handedness  difference,  and  our  phase  array  is 
defined  on  a  16x16  grid  whereas  the  Draper  phase  is  defined  on  a 
21x21  grid.  These  differences  made  it  impossible  to  show  that  the 
phase  estimation  algorithm  was  working  correctly. 


In  view  of  these  difficulties,  we  (Draper  and  EIKONIX  personnel) 
decided  to  attempt  a  reconciliation  of  the  differences  in  software 
and  to  have  EIKONIX  perform  another  round  of  testing,  this  time 
without  a  blindfold.  Thus  we  rewrote  our  software  to  produce  a  point 
spread  function  based  on  a  21x21  grid  for  the  phase  and  to  sort  out 
the  handedness  problem. 

We  buffered  the  21x21  grid  (in  the  spatial  frequency  domain) 
with  zeroes  to  form  a  64x64  grid.  Then  we  used  a  2-D,  64x64  FFT  to 
compute  a  coherent  impulse  response,  squared  the  samples,  and 
averaged  them  with  a  6. 5x6. 5  detector.  This  produced  a  9x9  array 
from  which  we  selected  an  8x8  array  to  be  used  as  the  observed  point 
spread  function. 

This  procedure  approximated  the  way  that  Draper  produced  the 
8x8  PSF  array.  The  major  difference  is  that  they  buffer  the  21x21 
array  to  128x128 ,  Fourier  transform  the  array,  square  the  samples, 
and  use  a  13x13  detector.  Also,  the  13x13  detector  size  is  realized 
by  frequency  domain  multiplication,  which  requires  another  pair  of 
128x128  Fourier  transforms .  Since  we  use  this  procedure  so  often 
during  our  phase  estimation  search,  we  could  not  afford  to  use  the 
larger  (128x128  vs.  64x64)  array  or  calculate  the  second  pair  of 
transforms. 
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Next,  we  sorted  out  the  handedness  problem.  This  was  done  by 
putting  the  actual  Draper  phase  into  our  algorithm  described  above. 

It  was  immediately  obvious  that  an  x-y  inversion  existed  and  that 
our  y  and  Draper's  y  were  defined  in  opposite  senses.  We  modified 
our  software  to  rectify  this  problem. 

In  Figure  73  we  show  the  Draper  phase  and  in  Figure  74  we  show 
the  Draper  PSF  and  the  EIKONIX  PSF.  The  differences  in  Figure  74 
are  so  slight  that  we  proceeded  with  confidence  into  the  phase 
estimation  algorithm. 

First  we  performed  phase  retrieval  on  the  noiseless  PSF  (the 
PSF  on  the  left  in  Figure  74) .  Since  we  are  no  longer  blind,  we  are 
able  to  see  how  well  the  estimation  is  done  at  each  step  in  the 
search.  The  results  are  summarized  in  Figure  75.  Starting  at  an 
all-zero  phase  assumption  we  find  the  best  6-parameter  phase  (an 
unimportant  constant  phase,  two  tilts,  and  3  quadratic  terms) ;  then, 
starting  at  this  solution,  we  find  the  best  10-parameter  phase,  etc., 
to  15  and  to  23  parameters. 

It  interesting  to  note  that  the  final  rms  phase  difference  between 
the  estimated  and  actual  phase  is  0.0542  waves.  In  discussion  with 
Draper  personnel,  we  find  that  this  is  approximately  the  same  error 
that  one  calculates  when  a  polynomial  of  this  size  is  fitted  directly 
to  the  phase  itself.  Thus,  we  conduce  that  phase  retrieval  on  the 
noiseless  PSF  was  successful.  Our  estimated  phase  is  shown  in  Figure 
76. 

Next  we  repeated  phase  retrieval  for  the  first  noisy  PSF, 

Figure  34(a).  This  was  derived  from  the  same  Draper  phase  but  noise 
is  added  to  the  PSF.  We  calculated  the  rms  phase  difference  between 
actual  and  estimated  by  normalizing  both  phases.  We  calculated  the 
constant  and  set  of  tilts  that  minimized  the  rms  phase  of  each,  then 
calculated  the  difference.  The  results  are  summarized  in  Figures  77 
and  78  . 
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Figure  77  shows  that  the  best  results  are  obtained  with  only 
10  parameters.  Apparently,  for  more  parameters  the  algorithm  starts 
to  use  the  extra  degrees  of  freedom  to  fit  the  noise.  Since  the 
metric  (MSE)  is  availalbe  for  inspection  and  since  the  expected 
noise  variance  may  be  known,  this  implies  that  one  could  continue 
using  more  parameters  in  the  search  until  the  expected  residual  MSE 
is  achieved.  There  the  search  would  be  terminated. 

We  note  that  for  this  noisy  PSF  the  expected  MSE  is 

MSE  =  (#  of  PSF  samples)  (noise  variance) 

=  64  (2%  of  DL-PSF  peak)2 
=  64  (0.02  *  835) 2 
=  17,848 

and  the  algorithm  yields  a  MSE  that  drops  from  18,277  to  15,494  as 
the  number  of  parameters  is  increased  from  6  to  10.  Thus  10  parameters 
is  a  reasonable  stopping  point. 

Finally  we  repeated  the  phase  retrieval  algorithm  on  the  second, 
noisy  PSF.  The  results  are  given  in  Figures  79  and  80. 

In  summary,  the  phase  retrieval  algorithm  was  reworked  so  that 
there  is  no  effective  miss-match  between  the  Draper  and  EIKONIX 
software.  The  polynomials  are  still  different  but  the  manner  in 
which  a  given  polynomial-generated  phase  is  converted  to  a  PSF  is 
the  same.  The  algorithm  then  yielded  a  phase  that  leaves  a  residual 
of  about  0.1  waves  between  the  actual  and  estimated  phases,  when 
the  input  PSF  is  undersampled  by  a  factor  of  4  in  both  x  and  y 
directions  and  when  there  is  a  2%  noise  on  the  measured,  undersampled 
PSF. 
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(a)  Noiseless  PSF  from  Draper  Labs 

(b)  EIKONIX  generated  PSF 
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Figure  75.  Phase  Retrieval  Summary  for  the  Noiseless  PSF 
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Figure  76.  Phase  Retrieval  Estimate  for  the  Noiseless  PSF 
(compare  to  Figure  73) 
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Figure  77.  (a)  First  noisy  PSF 

(b)  Estimated  PSF 

(c)  Phase  Retrieval  Summary 
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Figure  78.  Phase  Retrieval  Estimate  for  the  First  Noisy  PSF 
(compare  to  Figure  73) 
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Figure  79.  (a)  Input  PSF 

(b)  Estimated  Phase 

(c)  Phase  Retrieval  Summary 
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Figure  80.  Phase  Retrieval  Estimate  for  the  Second  Noisy  PSF 
(Compare  to  Figure  73) 
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